﻿# -*- coding: utf-8 -*-

# 任意路径线缆场线耦合效应之单线版__分段直线
import numpy as np
import environment_pool as env
def main():
	cos = np.cos
	sin = np.sin
	print("开始计算")
	eps = 2.220446049250313e-16
	# 常数
	c0 = 299792458
	mu0 = 4.0*np.pi*1e-7
	eps0 = 1.0/(mu0*c0*c0)
	j = 1.0j
	# 分析频率设置
	f = np.arange(env.Environment["起始频率"],env.Environment["终止频率"],env.Environment["步进频率"])#1e6:1e6:1000e6#频率
	flist = f.copy()
	Lenf = np.size(f)
	w = 2*np.pi*f
	beta = w/c0
	# 平面波参数[对标FEKO]
	E = env.Environment["E"]
	theta_p = env.Environment["theat"]
	psi_p = env.Environment["psi"]
	theta_E = env.Environment["eta"]
	
	expand_dims = lambda u:np.expand_dims(u, axis=-1)
	
	theta = (theta_p)*np.pi/180
	psi = (psi_p+90)*np.pi/180
	eta = (theta_E+180)*np.pi/180
	
	ex = cos(eta)*sin(theta)
	ey = sin(eta)*cos(psi)-cos(eta)*cos(theta)*sin(psi)
	ez = sin(eta)*sin(psi)+cos(eta)*cos(theta)*cos(psi)
	betax = beta*cos(theta)
	betay = beta*sin(theta)*sin(psi)
	betaz = beta*sin(theta)*cos(psi)
	# 向量维度变换
	betax = betax.reshape((1,1,Lenf))
	betay = betay.reshape((1,1,Lenf))
	betaz = betaz.reshape((1,1,Lenf))
	# 线缆路径导入
	Nt = int((env.Environment["线缆路径参数终止值"]-env.Environment["线缆路径参数起始值"])/env.Environment["线缆路径参数步长"])
	t = np.linspace(env.Environment["线缆路径参数起始值"], env.Environment["线缆路径参数终止值"], Nt) # 参数t的范围和点数
	if env.Environment["线缆路径文件"] is None:
		data = {
			"<线缆路径x公式>":env.Environment["线缆路径x公式"],
			"<线缆路径y公式>":env.Environment["线缆路径y公式"],
			"<线缆路径z公式>":env.Environment["线缆路径z公式"],
			"<线缆路径参数起始值>":env.Environment["线缆路径参数起始值"],
			"<线缆路径参数终止值>":env.Environment["线缆路径参数终止值"],
			"<Nt>":Nt,
		}
		print("<线缆路径x公式>__________"+env.Environment["线缆路径x公式"])
		from utils import Formula_cfg_and_run
		Formula_cfg_and_run(env.Environment["WorkPath"],data)
		result = np.loadtxt(env.Environment["WorkPath"]+'\\result.txt')
	else:
		result = np.loadtxt(env.Environment["线缆路径文件"],delimiter=',')
		Nt=np.size(result[:,0])	#利用线缆曲线点坐标拟合出线缆C表达式
		t=np.linspace(0,1,Nt)
	from spacecurve import SpaceCurve
	scurve = SpaceCurve(t)
	scurve.initial(result)
	fx,fy,fz = scurve.curve_fit()
	C = lambda t : np.vstack([fx(t) ,fy(t), fz(t)]) # [x(t) y(t) z(t)]
	
	# # 示例1：抛物线中心曲线
	# C = @(t) [-0.1 * t**2 zeros(size(t)) t] # [x(t) y(t) z(t)] x is the height
	p = 1
	h0 = 0
	L = 1
	# if env.Environment["线缆路径类型"]==0:
	#     C = lambda t : np.vstack([np.zeros(np.size(t)) ,np.zeros(np.size(t)), t]) # [x(t) y(t) z(t)]
	# elif env.Environment["线缆路径类型"]==1:
	#     C = lambda t : np.vstack([-p*(t-0.5*L)**2.0+0.25*p*L**2+h0 ,np.zeros(np.size(t)), t]) # [x(t) y(t) z(t)]
	# # C = @(t) [-p*(t-0.5*L)**2 zeros(size(t)) t] # [x(t) y(t) z(t)]
	
	# 示例2：圆形中心曲线（取消注释以下行）
	# C = @(t) [cos(t) sin(t) zeros(size(t))] # 半径为1的圆
	
	# 示例3：螺旋形中心曲线（取消注释以下行）
	# C = @(t) [0.5 * cos(t) 0.5 * sin(t) 0.2 * t]
	
	# 示例4：直线形中心曲线（取消注释以下行）
	# C = @(t) [zeros(size(t)) zeros(size(t)) t]
	
	# 计算Frenet标架（切向量T、法向量N、副法向量B）
	# 数值微分计算导数和二阶导数
	delta = 1e-5
	dC = lambda t :(C(t + delta) - C(t - delta)) / (2 * delta)
	
	# dC = @(t) diff(C(t))
	ddC = lambda t :(C(t + delta) - 2 * C(t) + C(t - delta)) / delta**2
	# ddC = @(t)diff(dC(t))
	
	# 初始化标架
	T = np.zeros((3, np.size(t)))
	N = np.zeros((3, np.size(t)))
	B = np.zeros((3, np.size(t)))
	
	
	for i in np.arange(np.size(t)):
		ti = t[i]
		dC_val = dC(ti)
		ddC_val = ddC(ti)
		
		# 切向量T
		T_val = dC_val / np.linalg.norm(dC_val)
		T_val = T_val.ravel()
		T[:, i] = T_val.ravel()
		# 副法向量B（叉积归一化）
		B_val = np.cross(dC_val.ravel(), ddC_val.ravel())
		B_val = (B_val)/(np.linalg.norm(B_val) + eps) # 避免除以零
		if np.linalg.norm(B_val)<1e-10:
			d = [0,0,1]
			B_val = np.cross(T_val, d)
			if np.linalg.norm(B_val)<1e-10:
				d = [0,1,0]
				B_val = np.cross(T_val, d)    
				B_val = B_val/np.linalg.norm(B_val)
			
				if np.linalg.norm(B_val)<1e-10:
					d = [1,0,0]
					B_val = np.cross(T_val, d)    
					B_val = B_val/np.linalg.norm(B_val)
				
					
		B[:, i] = B_val
		
		# 法向量N
		N_val = np.cross(B_val.ravel(), T_val.ravel())
		N[:, i] = N_val
	
	Temp = C(t)
	# 线缆参数设置
	Cable_Type = env.Environment["线缆类型"]
	a = env.Environment["内导体半径"]*1e-3# 单位mm
	b = env.Environment["线缆半径"]*1e-3# 单位mm
	epsr = env.Environment["介质层相对介电常数"]
	pos = np.array(env.Environment["中心起点坐标"],dtype=np.double).T*1e-3# 单位mm
	mur = 1
	s = env.Environment["线间距"]*1e-3# 单位mm
	pit = env.Environment["螺距"]*1e-3# 单位mm
	alpha = 1/np.sqrt((s/2)**2+(pit/(2*np.pi))**2)  # 角频率（控制缠绕密度）
	if Cable_Type == "双绞线" :
		Num = env.Environment["线缆数"]
		# 端接负载
	elif Cable_Type == "单线" :
		Num = env.Environment["线缆数"]
	elif Cable_Type == "四绞线" :
		Num = env.Environment["线缆数"]
		# 端接负载
	elif Cable_Type == "排线" :
		Num = env.Environment["线缆数"]
		# 端接负载
	R1 = env.Environment["左端端接负载"]*np.eye(Num)
	R2 = env.Environment["右端端接负载"]*np.eye(Num)
	
	x = np.zeros((Num,np.size(t)))
	y = np.zeros((Num,np.size(t)))
	z = np.zeros((Num,np.size(t)))
	for i in np.arange(Num):
		#P = Temp[i,:]+s/2*cos(alpha*t+i*np.pi)* N[i,:] + s/2*sin(alpha*t+i*np.pi)* B[i,:]
		if Cable_Type == "双绞线" :
			P = Temp+s/2*cos(alpha*t+i*np.pi)* N + s/2*sin(alpha*t+i*np.pi)* B+pos.reshape((-1,1))
		elif Cable_Type == "单线" :
			P = Temp+pos.reshape((-1,1))
		elif Cable_Type == "四绞线" :
			P = Temp+s/2*cos(alpha*t+i*np.pi/2)* N +s/2*np.sin(alpha*t+i*np.pi/2)* B+pos.reshape((-1,1))
		elif Cable_Type == "排线" :
			d_v = pos[:,0]-pos[:,i]
			P = Temp+d_v[0]*N+d_v[1]*B+pos[:,0].reshape((-1,1))
		x[i,:] = P[0,:]
		y[i,:] = P[1,:]
		z[i,:] = P[2,:]
	# print(x[:,:])
	# rw = 0.25e-3#线半径
	# 沿线场计算
	x0 = x
	y0 = y
	z0 = z
	x = expand_dims(x)
	y = expand_dims(y)
	z = expand_dims(z)
	Ex_Re = 2*E*ex*cos(betax*x)*cos(betay*y-betaz*z)# 矩阵维度解释：线缆编号*线段位置*频率
	Ex_Im = 2*E*ex*cos(betax*x)*sin(betay*y-betaz*z)
	Ey_Re = -2*E*ey*sin(betax*x)*sin(betay*y-betaz*z)
	Ey_Im = 2*E*ey*sin(betax*x)*cos(betay*y-betaz*z)
	Ez_Re = -2*E*ez*sin(betax*x)*sin(betay*y-betaz*z)
	Ez_Im = 2*E*ez*sin(betax*x)*cos(betay*y-betaz*z)
	x = x0
	y = y0
	z = z0
	# 设置初始值
	Num_Mesh = np.size(x[0,:])#线缆剖分段数
	Len_Slot = np.ones((Num,Num_Mesh-1))
	#是否含介质层判定
	epsr_diff = epsr-1
	# 等效线缆半径
	#epsr_diff==0
	re = (a**(1/epsr))*(b**(1-1/epsr))
	rm = (a**(mur))*(b**(1-mur))
	VsL = np.zeros((Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	VsR = np.zeros((Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	
	l = np.zeros((Num,Num,Num_Mesh-1))
	Zc = np.zeros((Num,Num,Num_Mesh-1),dtype=np.complex64)
	Zc_inv = np.zeros((Num,Num,Num_Mesh-1),dtype=np.complex64)
	p = np.zeros((Num,Num,Num_Mesh-1))
	c = np.zeros((Num,Num,Num_Mesh-1))
	phi11 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	phi12 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	phi21 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	phi22 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	
	for i in np.arange(Num_Mesh-1):
		Tmp = np.sqrt((x[:,i+1]-x[:,i])**2+(y[:,i+1]-y[:,i])**2+(z[:,i+1]-z[:,i])**2)
		Len_Slot[:,i] = Tmp.ravel()
		xi_x = (x[:,i+1]-x[:,i])/Len_Slot[:,i]# 单位方向向量的x分量
		xi_y = (y[:,i+1]-y[:,i])/Len_Slot[:,i]# 单位方向向量的y分量
		xi_z = (z[:,i+1]-z[:,i])/Len_Slot[:,i]# 单位方向向量的z分量
		tmp0 = np.outer(beta,Len_Slot[:,i])
		beta = beta.reshape((-1,1))
		
		Sigma_L = ((sin(tmp0)-tmp0)/(beta**2*sin(tmp0)))
		Sigma_R = ((sin(tmp0)-tmp0*cos(tmp0))/(beta**2*sin(tmp0)))
		Psi_L = ((cos(tmp0)-np.ones((Lenf,Num)))/(beta*sin(tmp0)))
		Psi_R = -((cos(tmp0)-np.ones((Lenf,Num)))/(beta*sin(tmp0)))
		ax = (Ex_Re[:,i+1,:]-Ex_Re[:,i,:])/expand_dims(Len_Slot[:,i])
		ax = ax.reshape((Num,Lenf))
		ay = (Ey_Re[:,i+1,:]-Ey_Re[:,i,:])/expand_dims(Len_Slot[:,i])
		ay = ay.reshape((Num,Lenf))
		az = (Ez_Re[:,i+1,:]-Ez_Re[:,i,:])/expand_dims(Len_Slot[:,i])
		az = az.reshape((Num,Lenf))
		bx = Ex_Re[:,i,:]
		bx = bx.reshape((Num,Lenf))
		by = Ey_Re[:,i,:]
		by = by.reshape((Num,Lenf))
		bz = Ez_Re[:,i,:]
		bz = bz.reshape((Num,Lenf))
		cx = (Ex_Im[:,i+1,:]-Ex_Im[:,i,:])/expand_dims(Len_Slot[:,i])
		cx = cx.reshape((Num,Lenf))
		cy = (Ey_Im[:,i+1,:]-Ey_Im[:,i,:])/expand_dims(Len_Slot[:,i])
		cy = cy.reshape((Num,Lenf))
		cz = (Ez_Im[:,i+1,:]-Ez_Im[:,i,:])/expand_dims(Len_Slot[:,i])
		cz = cz.reshape((Num,Lenf))
		dx = Ex_Im[:,i,:]
		dx = dx.reshape((Num,Lenf))
		dy = Ey_Im[:,i,:]
		dy = dy.reshape((Num,Lenf))
		dz = Ez_Im[:,i,:]
		dz = dz.reshape((Num,Lenf))
		xi_x = expand_dims(xi_x)
		xi_y = expand_dims(xi_y)
		xi_z = expand_dims(xi_z)
		# 等效耦合源
		VsL_l = xi_x*((ax+j*cx)*Sigma_L.T+(bx+j*dx)*Psi_L.T)+xi_y*((ay+j*cy)*Sigma_L.T+(by+j*dy)*Psi_L.T)+xi_z*((az+j*cz)*Sigma_L.T+(bz+j*dz)*Psi_L.T)
		VsR_l = xi_x*((ax+j*cx)*Sigma_R.T+(bx+j*dx)*Psi_R.T)+xi_y*((ay+j*cy)*Sigma_R.T+(by+j*dy)*Psi_R.T)+xi_z*((az+j*cz)*Sigma_R.T+(bz+j*dz)*Psi_R.T)
		VsL_v = -(Ex_Re[:,i,:]+j*Ex_Im[:,i,:])*expand_dims(x[:,i])
		VsL_v = VsL_v.reshape((Num,Lenf))
		VsR_v = -(Ex_Re[:,i+1,:]+j*Ex_Im[:,i+1,:])*expand_dims(x[:,i+1])
		VsR_v = VsR_v.reshape((Num,Lenf))
		VsL[:,:,i] = VsL_l+VsL_v
		VsR[:,:,i] = VsR_l+VsR_v# (Num,freq,Mesh)
		# 单位长度参数
		for ii in np.arange(Num) :
			if Cable_Type == "四绞线" :
				if ii == 0:
					l[ii,ii,i] = (mu0/(2*np.pi))*(np.log(2*(Temp[0,i]+pos[0])/rm)-s**2/(16*(Temp[0,i]+pos[0])**2))
				else:
					l[0,ii,i] = (mu0/(2*np.pi))*(np.log(2*(Temp[0,i]+pos[0])/(s*np.abs(np.sin(ii*np.pi/4))))-s**2/(16*(Temp[0,i]+pos[0])**2)*np.cos(ii*np.pi/2))
			else:
				for jj in np.arange(Num) :
					if ii == jj:
						if Cable_Type == "双绞线" :
							l[ii,jj,i] = (mu0/(2*np.pi))*np.log(2*Temp[0,i]+pos[0]/rm)
							p[ii,jj,i] = (1/(2*np.pi*eps0))*np.log(2*Temp[0,i]+pos[0]/re)
						elif Cable_Type == "单线" or Cable_Type == "排线" :
							l[ii,jj,i] = (mu0/(2*np.pi))*np.log(2*x[ii,i]/rm)
							p[ii,jj,i] = (1/(2*np.pi*eps0))*np.log(2*x[ii,i]/re)
					else:
						if Cable_Type == "双绞线" :
							l[ii,jj,i] = (mu0/(2*np.pi))*(np.log(2*(Temp[0,i]+pos[0])/s)+s**2/(16*(Temp[0,i]+pos[0])**2))
							p[ii,jj,i] = (1/(2*np.pi*eps0))*(np.log(2*(Temp[0,i]+pos[0])/s)+s**2/(16*(Temp[0,i]+pos[0])**2))
						elif Cable_Type == "单线" or Cable_Type == "排线" :
							l[ii,jj,i] = (mu0/(4*np.pi))*np.log(1+4*x[ii,i]*x[jj,i]/((x[ii,i]-x[jj,i])**2+(y[ii,i]-y[jj,i])**2+(z[ii,i]-z[jj,i])**2))
							p[ii,jj,i] = (1/(4*np.pi*eps0))*np.log(1+4*x[ii,i]*x[jj,i]/((x[ii,i]-x[jj,i])**2+(y[ii,i]-y[jj,i])**2+(z[ii,i]-z[jj,i])**2))
		if epsr_diff>0:
			c[:,:,i]= np.linalg.pinv(p[:,:,i])
		else:
			if Cable_Type == "四绞线" :
				l_ind = np.array([[l[0,0,i],l[0,1,i],l[0,2,i],l[0,3,i]],
					[l[0,1,i],l[0,0,i],l[0,1,i],l[0,2,i]],
					[l[0,2,i],l[0,1,i],l[0,0,i],l[0,1,i]],
					[l[0,1,i],l[0,2,i],l[0,1,i],l[0,0,i]]])
				Zc[:,:,i] = c0*l_ind
			else:
				Zc[:,:,i] = c0*l[:,:,i]
			Zc_inv[:,:,i] = np.linalg.pinv(Zc[:,:,i])
	for i in np.arange(Lenf):
		for ii in np.arange(Num_Mesh-1):
			if(epsr_diff>0):
				from scipy.linalg import eig
				Z = j*w[i]*l[:,:,ii]
				Y = j*w[i]*c[:,:,ii]
				Y_inv = np.linalg.pinv(Y)
				Gam,TI=eig(Y@Z)
				TI_inv = np.linalg.pinv(TI)
				Gam = np.sqrt(Gam*np.eye(Num))
				Zc[:,:,ii] = Y_inv@TI@Gam@TI_inv
				Zc_inv[:,:,ii] = np.linalg.pinv(Zc[:,:,ii])
			# 传输链参数矩阵
			tmp0 = beta[i]*Len_Slot[0,ii]
			tmp0 = tmp0[0]
			phi11[:,:,i,ii] = cos(tmp0)*np.eye(Num,dtype=np.float64)
			phi12[:,:,i,ii] = -1j*sin(tmp0)*Zc[:,:,ii]
			phi21[:,:,i,ii] = -1j*sin(tmp0)*Zc_inv[:,:,ii]
			phi22[:,:,i,ii] = cos(tmp0)*np.eye(Num,dtype=np.float64)
	def write_curve3D(fileName,data):
		ndata = np.size(data[:,0])
		f = open(fileName,"w")
		f.write(f"{ndata}\n")
		np.savetxt(f,data)
		f.close()
	write_curve3D(f'{env.Environment["中心线缆名称"]}.txt',Temp.T)
	for i in np.arange(Num):
		data = np.zeros((np.size(t),3))
		data[:,0] = x[i,:]
		data[:,1] = y[i,:]
		data[:,2] = z[i,:]
		write_curve3D(f'{env.Environment["线缆集名称"]}{i}.txt',data)
	Phi11 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	Phi11[:,:,:,0] =phi11[:,:,:,Num_Mesh-2]
	Phi12 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	Phi12[:,:,:,0] =phi12[:,:,:,Num_Mesh-2]
	Phi21 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	Phi21[:,:,:,0] = phi21[:,:,:,Num_Mesh-2]
	Phi22 = np.zeros((Num,Num,Lenf,Num_Mesh-1),dtype=np.complex128)
	Phi22[:,:,:,0] =phi22[:,:,:,Num_Mesh-2]
	Cons_V = np.zeros((Num,Lenf),dtype=np.complex128)
	Cons_I = np.zeros((Num,Lenf),dtype=np.complex128)
	np.set_printoptions(precision=20)
	for i in np.arange(Lenf):
		Cons_V[:,i] = VsR[:,i,Num_Mesh-2]+Phi11[:,:,i,0]@(-VsL[:,i,Num_Mesh-2]+VsR[:,i,Num_Mesh-3])
		Cons_I[:,i] = Phi21[:,:,i,0]@(-VsL[:,i,Num_Mesh-2]+VsR[:,i,Num_Mesh-3])
		for ii in np.arange(1,Num_Mesh-2):
			Phi11[:,:,i,ii] = Phi11[:,:,i,ii-1]@phi11[:,:,i,Num_Mesh-ii-2]+Phi12[:,:,i,ii-1]@phi21[:,:,i,Num_Mesh-ii-2]
			Phi12[:,:,i,ii] = Phi11[:,:,i,ii-1]@phi12[:,:,i,Num_Mesh-ii-2]+Phi12[:,:,i,ii-1]@phi22[:,:,i,Num_Mesh-ii-2]
			Phi21[:,:,i,ii] = Phi21[:,:,i,ii-1]@phi11[:,:,i,Num_Mesh-ii-2]+Phi22[:,:,i,ii-1]@phi21[:,:,i,Num_Mesh-ii-2]
			Phi22[:,:,i,ii] = Phi21[:,:,i,ii-1]@phi12[:,:,i,Num_Mesh-ii-2]+Phi22[:,:,i,ii-1]@phi22[:,:,i,Num_Mesh-ii-2]
			Cons_V[:,i] = Cons_V[:,i]+Phi11[:,:,i,ii]@(-VsL[:,i,Num_Mesh-ii-2]+VsR[:,i,Num_Mesh-ii-3])
			Cons_I[:,i] = Cons_I[:,i]+Phi21[:,:,i,ii]@(-VsL[:,i,Num_Mesh-ii-2]+VsR[:,i,Num_Mesh-ii-3])
	I0 = np.zeros((Num,Lenf),dtype=np.complex64)
	V0 = np.zeros((Num,Lenf),dtype=np.complex64)
	IL = np.zeros((Num,Lenf),dtype=np.complex64)
	VL = np.zeros((Num,Lenf),dtype=np.complex64)
	for i in np.arange(Lenf):
		Phi11[:,:,i,Num_Mesh-2] = Phi11[:,:,i,Num_Mesh-3]@phi11[:,:,i,0]+Phi12[:,:,i,Num_Mesh-3]@phi21[:,:,i,0]
		Phi12[:,:,i,Num_Mesh-2] = Phi11[:,:,i,Num_Mesh-3]@phi12[:,:,i,0]+Phi12[:,:,i,Num_Mesh-3]@phi22[:,:,i,0]
		Phi21[:,:,i,Num_Mesh-2] = Phi21[:,:,i,Num_Mesh-3]@phi11[:,:,i,0]+Phi22[:,:,i,Num_Mesh-3]@phi21[:,:,i,0]
		Phi22[:,:,i,Num_Mesh-2] = Phi21[:,:,i,Num_Mesh-3]@phi12[:,:,i,0]+Phi22[:,:,i,Num_Mesh-3]@phi22[:,:,i,0]
		Cons_V[:,i] = Cons_V[:,i]+Phi11[:,:,i,Num_Mesh-2]@(-VsL[:,i,0])
		Cons_I[:,i] = Cons_I[:,i]+Phi21[:,:,i,Num_Mesh-2]@(-VsL[:,i,0])
		# 终端电压电流
		I0[:,i] = np.linalg.solve(Phi12[:,:,i,Num_Mesh-2]-Phi11[:,:,i,Num_Mesh-2]@R1+R2@Phi21[:,:,i,Num_Mesh-2]@R1-R2@Phi22[:,:,i,Num_Mesh-2],R2@Cons_I[:,i]-Cons_V[:,i])
		V0[:,i] = -R1@I0[:,i]
		IL[:,i] = Cons_I[:,i]+Phi21[:,:,i,Num_Mesh-2]@V0[:,i]+Phi22[:,:,i,Num_Mesh-2]@I0[:,i]
		VL[:,i] = R2@IL[:,i]
	
	# I0_dBuA = 120+20*log10(abs(I0))
	# IL_dBuA = 120+20*log10(abs(IL))
	# V0_dBuA = 120+20*log10(abs(V0))
	# VL_dBuA = 120+20*log10(abs(VL))
	V0_dBV = 20*np.log10(np.abs(V0))
	VL_dBV = 20*np.log10(np.abs(VL))
	I0_dBV = 20*np.log10(np.abs(I0))
	IL_dBV = 20*np.log10(np.abs(IL))
	
	for ii in np.arange(Num):
		data = np.zeros((Lenf,13))
		data[:,0] = flist
		data[:,1] = V0_dBV[ii,:]
		data[:,2] = VL_dBV[ii,:]
		data[:,3] = I0_dBV[ii,:]
		data[:,4] = IL_dBV[ii,:]
	
		data[:,5] = np.angle(V0[ii,:]) 
		data[:,6] = np.angle(VL[ii,:]) 
		data[:,7] = np.angle(I0[ii,:]) 
		data[:,8] = np.angle(IL[ii,:]) 
	
		data[:,9] =  np.abs(V0[ii,:]) 
		data[:,10] = np.abs(VL[ii,:]) 
		data[:,11] = np.abs(I0[ii,:]) 
		data[:,12] = np.abs(IL[ii,:]) 
		
		f = open(f'{env.Environment["结果集名称"]}{ii}.txt',"w")
		f.write("Freq V0_dBV VL_dBV I0_dBV IL_dBV V0_phase VL_phase I0_phase IL_phase V0_abs VL_abs I0_abs IL_abs\n")
		np.savetxt(f,data)
		f.close()
if __name__ == "__main__":
    # 作为脚本直接运行时调用 main()
    main()
